/*
 *  MWSolve.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 10-07-12.
 *  Copyright 2010 MikBarr Studios. All rights reserved.
 *
 */

#include "MW.h"
#include "RKModMaxwell.h"
#include <cmath>
#include <sys/time.h>
#include <iostream>
#include <fstream>
using namespace T3;

Scalar adiff(Scalar,real);

void MW::Solve() {
	std::ofstream info((id + "/info.txt").c_str());
	Params(info);
	info.close();
	
	Badge->PrintStatus(id);
    
	Badge->PrintStatus("Forming solution.");
	real dt = s*O.d();
	int N = int(T/dt);
	real t=0;
	real I0 = I;
    
	long long int StoSize = 12LL*8*N*(7*n*n+1);
    
	Badge->PrintStatus((string)"Estimated storage usage: " + bytes(StoSize));
    
#define Phi M.S(0)
#define Psi M.S(1)
#define E M.V(0)
#define B M.V(1)
    
	Set M(O,3,2,2);
	Phi.id = "Phi";
	Psi.id = "Psi";
	E.id = "E";
	E(0).id = "E0";
	E(1).id = "E1";
	E(2).id = "E2";
	B.id = "B";
	B(0).id = "B0";
	B(1).id = "B1";
	B(2).id = "B2";
    
	Badge->PrintStatus("Initial data");
	/* Plane wave in z */
    /*	Scalar R = smoother(x(2), 1.);
     Scalar dR = diff(R)(2);
     FORO {
     E(0)[o] = omega*(R[o]*cos(omega*x(2)[o])-dR[o]*sin(omega*x(2)[o]));
     E(1)[o] = omega*(-R[o]*sin(omega*x(2)[o])-dR[o]*cos(omega*x(2)[o]));
     B(0)[o] = -E(1)[o];
     B(1)[o] = E(0)[o];
     }
     if ( c != 0 ) {
     Badge->PrintStatus("Initializing violation");
     FORO B(0)[o] += c*R[o];
     }*/
    
	/* Spherical wave */
    /*	real rpos = 1.;
     real rdel = 5*O.d();
     Scalar rmin(O), rmax(O), rinv(O);
     FORO rmin[o] = x.rad(o)-(rpos-rdel/2);
     FORO rmax[o] = x.rad(o)-(rpos+rdel/2);
     FORO rinv[o] = ( x.rad(o) == 0 ) ? 0. : 1./x.rad(o);
     Scalar R(O), dR(O), d2R(O);
     FORO R[o] = smoother(rmin[o], rdel)*smoother(rmax[o], rdel);
     FORO dR[o] = dsmoother(rmin[o], rdel)*smoother(rmax[o], rdel)+smoother(rmin[o], rdel)*dsmoother(rmax[o], rdel);
     FORO d2R[o] = d2smoother(rmin[o], rdel)*smoother(rmax[o], rdel)+2*dsmoother(rmin[o], rdel)*dsmoother(rmax[o], rdel)+smoother(rmin[o], rdel)*d2smoother(rmax[o], rdel);
     
     Scalar F1(O), F2(O), F3(O);
     FORO F1[o] = d2R[o]*pow(rinv[o],2)+dR[o]*pow(rinv[o],3);
     FORO F2[o] = d2R[o]*pow(rinv[o],3)+3*dR[o]*pow(rinv[o],4)+3*R[o]*pow(rinv[o],5);
     FORO F3[o] = d2R[o]*rinv[o]+dR[o]*pow(rinv[o],2)+R[o]*pow(rinv[o],3);
     
     FORO E(1)[o] = F1[o]*x(2)[o];
     FORO E(2)[o] = -F1[o]*x(1)[o];
     FORa FORO B(a)[o] = x(0)[o]*F2[o]*x(a)[o];
     FORO B(0)[o] -= F3[o];
     
     real vio = c*pow(10.,-3)*max(F1);
     FORO E(0)[o] = vio*R[o]; */
    
    /* Neutron star */
    real mu = 8000; // Magnetic field strength
	real sint = pole;
	real cost = 1-pole;
    // e_x = -e_phi, e_y = e_theta, e_z = e_r
    /* FIXME Big-box */
    /*real R = 1; // Radius of star
    PFOR(j,O(1)) FOR(i,O(0)) {
        real s = sqrt(R*R-x(1)[o]*x(1)[o]-x(0)[o]*x(0)[o]);
        B(1)(i,j,0) = -mu*sin(theta)/(s*s*s);
        B(2)(i,j,0) = 2*mu*cos(theta)/(s*s*s);
        E(1)(i,j,0) = -2*Omega*mu*sin(theta)*cos(theta)/(s*s);
        E(2)(i,j,0) = Omega*mu*pow(sin(theta),2.)/(s*s);
    }*/
    /* Small-box */
    real R = x(2)(0,0,0);
    FOR(i,O(0)) FOR(j,O(1)) FOR(k,O(2)) {
        real r = x(2)(i,j,k);
        B(1)(i,j,k) = mu*sint/(r*r*r);
        B(2)(i,j,k) = 2*mu*cost/(r*r*r);
        E(1)(i,j,k) = -2*Omega*mu*sint*cost/(r*r);
        E(2)(i,j,k) = Omega*mu*sint*sint/(r*r);
    }
    PFORO FOR(a,3) E(a)[o] += 0.2*B(a)[o];
    
	RKModMaxwell RKevo;
	RKevo.K = K;
	RKevo.I = I;
	RKevo.epsilon = epsilon;
    RKevo.R = R;
    RKevo.mu = mu;
    RKevo.pole = pole;
    RKevo.Omega = Omega;
    
	Badge->PrintStatus("Evolution loop");
	struct timeval start = gettime();
	for ( int n=0; t<T+dt/2 and !stopit; n++ ) {
/*		if ( fabs(t-floor(t)) < dt and I != 0 )
			RKevo.I = I0/(l2(E*E+B*B)+1);*/
		/*if ( n == N/10 )
			PFORO FOR(a,3) E(a)[o] += O.d()O.d()*B(a)[o];*/
		struct timeval first = gettime();
		Print(t, M);
		Badge->PrintTime(t);
		t += dt;
		RKevo.Evolve(M, dt);
		Badge->PrintStatus((string)"Time remaining: " + timer(first,N-n) + " at coordinate " + ftos(t),false);
		Badge->Progress(t/T);
	}
	Badge->PrintStatus((string)"\nTotal elapsed time: " + timer(start));
}

Scalar adiff(Scalar f, real F0) {
	Scalar F(f.O);
	FOR(i,F.O(0)) FOR(j,F.O(1)) {
		F(i,j,0) = F0;
		FOR(k,F.O(2)-1)
        F(i,j,k+1) = F(i,j,k) + F.O.d()*(f(i,j,k)+f(i,j,k+1))/2;
	}
	return F;
}

